M2.855 · Modelos avanzados de minería de datos · PEC2

2021-1 · Máster universitario en Ciencia de datos (Data science)

Estudios de Informática, Multimedia y Telecomunicación

 

PEC 2: Métodos no supervisados

A lo largo de esta práctica veremos como aplicar distintas técnicas no supervisadas así como algunas de sus aplicaciones reales:

• Clustering con distintas estrategias: k-means y regla del codo, basadas en densidad y jerárquicas.
• Aplicación: generación de imágenes con reducción de dimensionalidad. PCA y UMAP.
• Aplicación: identificación de puntos de interés turísticos.

Nombre y apellidos: Mario Ubierna San Mamés

---

Para ello vamos a necesitar las siguientes librerías:

import random

import umap
import numpy as np
import pandas as pd
import sklearn
from sklearn import cluster        # Algoritmos de clustering.
from sklearn import datasets       # Crear datasets.
from sklearn import decomposition  # Algoritmos de reduccion de dimensionalidad.

# Visualizacion.
import matplotlib
import matplotlib.pyplot as plt

%matplotlib inline

1. Métodos de clustering (4 puntos)

Este ejercicio trata de explorar distintas técnicas de agrupamiento ajustándolas a distintos conjuntos de datos.

El objetivo es doble: entender la influencia de los parámetros en su comportamiento, y conocer sus limitaciones en la búsqueda de estructuras de datos.

Generación de los conjuntos de datos

X_blobs, y_blobs = datasets.make_blobs(n_samples=1000, n_features=2, centers=4, cluster_std=1.6, random_state=42)
X_moons, y_moons = datasets.make_moons(n_samples=1000, noise=.07, random_state=42)
X_circles, y_circles = datasets.make_circles(n_samples=1000, factor=.5, noise=.05, random_state=42)

Cada dataset tiene 2 variables: una variable X que contiene 2 features (columnas) y tantas filas como muestras. Y una variable y que alberga las etiquetas que identifican cada cluster.

A lo largo del ejercicio no se usará la variable y (sólo con el objetivo de visualizar). El objetivo es a través de los distintos modelos de clustering conseguir encontrar las estructuras descritas por las variables y.

# NO ME FUNCIONA EN LOCAL (muere el kernel), SÍ EN GOOGLE COLAB
fig, axis = plt.subplots(2, 3, figsize=(13, 8))
for i, (X, y, ax, name) in enumerate(zip([X_blobs, X_moons, X_circles] * 2,
                                         [None] * 3 + [y_blobs, y_moons, y_circles],
                                         axis.reshape(-1),
                                         ['Blob', 'Moons', 'Circles'] * 2)):
    ax.set_title('Dataset: {}, '.format(name) + ('lo que analizarás' if i < 3 else 'los grupos a encontrar'))
    ax.scatter(X[:,0], X[:,1], s=15, c=y, alpha=.3, cmap='jet')
    ax.axis('equal')
plt.tight_layout()

1 a. K-means

En este apartado se pide probar el algoritmo k-means sobre los tres datasets presentados anteriormente ajustando con los parámetros adecuados y analizar sus resultados.

X, y = X_blobs, y_blobs

Para estimar el número de clusters a detectar por k-means. Una técnica para estimar $k$ es, como se explica en la teoría:

Los criterios anteriores (minimización de distancias intra grupo o maximización de distancias inter grupo) pueden usarse para establecer un valor adecuado para el parámetro k. Valores k para los que ya no se consiguen mejoras significativas en la homogeneidad interna de los segmentos o la heterogeneidad entre segmentos distintos, deberían descartarse.

Lo que popularmente se conocer como regla del codo.

Primero es necesario calcular la suma de los errores cuadráticos (SSE) que consiste en la suma de todos los errores (distancia de cada punto a su centroide asignado) al cuadrado.

$$SSE = \sum_{i=1}^{K} \sum_{x \in C_i} euclidean(x, c_i)^2$$

Donde $K$ es el número de clusters a buscar por k-means, $x \in C_i$ son los puntos que pertenecen a i-ésimo cluster, $c_i$ es el centroide del cluster $C_i$ (al que pertenece el punto $x$), y $euclidean$ es la distancia euclídea.

Este procedimiento realizado para cada posible valor $k$, resulta en una función monótona decreciente, donde el eje $x$ representa los distintos valores de $k$, y el eje $y$ el $SSE$. Intuitivamente se podrá observar un significativo descenso del error, que indicará el valor idóneo de $k$.

Se pide realizar la representación gráfica de la regla del codo junto a su interpretación, utilizando la librería matplotlib y la implementación en scikit-learn de k-means.

Implementación: cálculo y visualización de la regla del codo en el dataset Blobs.

def sse_kmeans(X, k): 
    sse = []
    for i in range(1, k):
        kmeans = cluster.KMeans(n_clusters=i)
        kmeans.fit(X)
        sse.append(kmeans.inertia_)
    return sse

def elbow_kmeans(sse, k):
    fig, ax = plt.subplots(figsize=(13, 8))
    ax.plot(range(1, k), sse)
    plt.title("Estimación del número de clústers")
    plt.xlabel("Número de clústers (k)")
    plt.ylabel("SSE")
    plt.show()

# Cálculo
k = 11
sse = sse_kmeans(X, k)

# Visualación
elbow_kmeans(sse, k)

C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4.
  warnings.warn(

Análisis: ¿Qué se interpreta en la gráfica? ¿Cómo podría mejorarse la elección de $k$?.

De la gráfica anterior podemos interpretar lo siguiente:

• El número de clústers óptimo para este dataset es de 4, es decir, elegimos como número de clústers (k) donde se produce el "codo" en la función, esto quiere decir donde la suma de los errores cuadráticos se atenúa (donde la distancia de cada punto a su centroide es menor).

Cabe destacar, que se podría elegir 3 clústers, pero este valor no sería el óptimo ya que la suma de los errores cuadráticos es mayor.

El SSE (suma de los errores cuadráticos), calcula la distancia de cada punto a su centroide, por lo que busca minimizar dicha distancia, es decir, que los clústers sean lo más compactos posibles. Otra forma de mejorar la elección de k es medir la distancia entre centroides, cuanto más distancia haya entre centroides más diferentes van a ser los clústers.

Por lo tanto, se podría hacer el cálculo de maximizar la distancia entre centroides para también mejorar la elección de k, en este caso porque en la anterior ejecución queda claro que lo mejor son 4 clústers, pero hacer este cálculo podría ayudarnos a si es mejor 3 o 5 clústers por ejemplo.

Implementación: cálculo y visualización de los grupos en el dataset Blobs.

def clusters_view(X, labels):
    fig, ax = plt.subplots(figsize=(13, 8))
    plt.title("Clústers detectados")
    
    colores = ["blue", "red", "green", "orange", "black", "grey"]
    labels_color = []
    for i in labels:
        labels_color.append(colores[i])
    ax.scatter(X[:,0], X[:,1], c=labels_color, alpha=.3)

# Cálculo
modelo_kmeans_blobs = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_blobs.fit(X)

# Visualización
clusters_view(X, modelo_kmeans_blobs.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

Como resultado de la anterior ejecución, podemos ver que tenemos 4 grupos tal y como era de esperar, ya que al inicializar kmeans el valor de k era 4.

Mencionar que los grupos obtenidos son correctos, es decir, tenemos 4 colores (4 grupos) y los puntos de un mismo color están más próximos al centro del clúster del mismo color que al centro de los otros grupos.

El resultado al ejecutar kmeans nos generá los siguientes grupos:

• Grupo arriba a la izquierda: se caracteriza porque el valor de la primera columna del dataset es pequeño, y el valor de la segunda columna es elevado.
• Grupo arriba a la derecha: se caracteriza porque el valor de la primera columna es medio y el de la segunda columna es elevado.
• Grupo abajo a la izquierda: tanto la primera columna como la segunda tienen valores pequeños.
• Grupo abajo a la derecha: en la primera columna tienen un elevado valor, mientras que en la segunda tienen un valor medio

X, y = X_moons, y_moons

Implementación: cálculo y visualización de la regla del codo en el dataset Moons.

# Cálculo
k = 11
sse = sse_kmeans(X, k)

# Visualación
elbow_kmeans(sse, k)

C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4.
  warnings.warn(

Análisis: ¿Qué se interpreta en la gráfica? ¿Cómo podría mejorarse la elección de $k$?.

Según la interpretación de la gráfica el valor óptimo del número de clústers podría ser 4 o 5, en nuestro caso vamos a elegir 4 ya que éste minimiza el error cuadrático de los residuos y generaliza algo más, por lo que a priori esta es la mejor opción.

Al igual que sucedía en el ejercicio anterior, podríamos mejorar la elección de k si calculásemos también la distancia entre los centroides, de tal forma que fuera la máxima, así se conseguría tener clústers bien diferentes unos de los otros a parte de ser homogéneos.

También se podría mejorar la elección de k probando diferentes modelos con diferentes k, en este caso no se sabe a priori si un modelo con 4 grupos es mejor que 5, por lo que se puede generar diferentes modelos y ver los resultado que proporciona.

Implementación: cálculo y visualización de los grupos en el dataset Moons.

# Cálculo
modelo_kmeans_moons = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_moons.fit(X)

# Visualización
clusters_view(X, modelo_kmeans_moons.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

En este caso, vemos que el modelo generado no es bueno, es decir, no consigue clasificar bien las instancias del dataset. Esto se debe principalmente a que el dataset es muy variado.

Aún así vemos que realmente sí que hay grupos, para ser más exactos hay dos, dos medias lunas, pero kmeans no es capaz de detectar bien los grupos básicamente porque el problema ha cambiado, los centroides de cada clúster no están claros, aunque sabemos que hay dos clústers, en resumen es un problema de densidad no de distancia. Kmeans no es capaz de elegir dos buenos centroides porque presupone que los clústers tienen forma convexa, no funciona bien si los clústers son no convexos, es decir, calcula los clústers a partir de la distancia no de la densidad.

X, y = X_circles, y_circles

Implementación: cálculo y visualización de la regla del codo en el dataset Circles.

# Cálculo
k = 11
sse = sse_kmeans(X, k)

# Visualación
elbow_kmeans(sse, k)

C:\Users\mario\anaconda3\envs\uoc20211pec1\lib\site-packages\sklearn\cluster\_kmeans.py:881: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=4.
  warnings.warn(

Análisis: ¿Qué se interpreta en la gráfica? ¿Cómo podría mejorarse la elección de $k$?.

De la gráfica anterior podemos concluir que el mejor número de clústers puede ser 4 o 5, ya que es en estos puntos donde se reduce más el error cuadrático de los residuos y así vez conseguimos generalizar. En nuestro caso, se ha decidido elegir un valor de k igual a 4, ya que conseguimos reducir el error y generalizar más.

Al igual que en los apartados anteriores, podemos mejorar la elección de k a partir de maximizar la distancia entre los centroides de cada grupo, o también ejecutando varios modelos y ver cuál nos produce un resultado mejor.

Implementación: cálculo y visualización de los grupos en el dataset Circles.

# Cálculo
modelo_kmeans_circles = cluster.KMeans(n_clusters=4, random_state=0)
modelo_kmeans_circles.fit(X)

# Visualización
clusters_view(X, modelo_kmeans_circles.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

En este caso pasa lo mismo que en el ejericio anterior, kmeans calcula los clústers a partir de la distancia de cada punto al centroide y no por densidad, por lo que no es capaz de distinguir bien cuántos clústers necesita y cuáles son éstos.

1 b. Algoritmos basados en densidad: DBScan

En este apartado se pide aplicar clustering por densidad como DBSCAN a los datasets anteriores para detectar los grupos subyacentes.

X, y = X_blobs, y_blobs

Implementación: prueba la implementación de DBSCAN en scikit-learn jugando con los parámetros eps y min_samples para encontrar los grupos (y outliers) del dataset Blobs.

# http://exponentis.es/ejemplo-de-uso-de-dbscan-en-python-para-deteccion-de-outliers
# Método que devuelve la regla del codo para un epsilón adecuado a nuestro problema
def elbow_dbscan(X, min_samples):
    estimator = decomposition.PCA(n_components = 2)
    X_pca = estimator.fit_transform(X)
    dist = sklearn.neighbors.DistanceMetric.get_metric('euclidean')
    matsim = dist.pairwise(X_pca)
    minPts  = min_samples # Fijamos el parámetro minPts
    A = sklearn.neighbors.kneighbors_graph(X_pca, minPts, include_self=False)
    Ar = A.toarray()
    seq = []
    for i,s in enumerate(X_pca):
        for j in range(len(X_pca)):
            if Ar[i][j] != 0:
                seq.append(matsim[i][j])
    seq.sort()
    plt.subplots(figsize=(13, 8))
    plt.title("Estimación de epsilon")
    plt.xlabel("Puntos")
    plt.ylabel("Distancia")
    plt.plot(seq)
    plt.show()

Se ha dictaminado que el min_samples sea 15 ya que era lo que mejor funcionaba con nuestro modeloo, es decir, identificaba bien los 4 clústers, ya que de lo contrario identificaba solo 3. A partir de min_samples se calcula con la función anterior el epsilon estimado.

# Epsilón adecuado
elbow_dbscan(X, 15)

# Cálculo
modelo_dbscan_blobs = cluster.DBSCAN(eps=1, min_samples=15)
modelo_dbscan_blobs.fit(X)

# Visualización
clusters_view(X, modelo_dbscan_blobs.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

En este caso vemos que sí que detecta bien los clústers al igual que antes, con la diferencia de que ahora los clústers se calculan a partir de la densidad y no de la distancia. Con eps indicamos el radio de vecindad de un punto, y con min_samples definimos la cantidad mínima de vecinos que hay dentro del radio de vecindad para considerar el punto dentro del clúster.

Cabe destacar que ahora no todos los puntos pertenecen a un clúster en sí, es decir, hay puntos grises que no forman parte de ningún clúster, en otras palabras tenemos outliers. Estos puntos se caracterizan porque no cumplen tanto la distancia como el número mínimo de vecinos en dicha distancia, por lo tanto, ahora a diferencia de kmeans podemos ver que no todos los puntos deberían pertenecer a un clúster, ya que no presentan las mismas similitudes.

X, y = X_moons, y_moons

Implementación: prueba la implementación de DBScan jugando con los parámetros eps y min_samples para encontrar los grupos (y outliers) del dataset Moons.

En este caso se ha dictaminado que min_samples sea 5 ya que era lo que mejor funcionaba para nuestro modelo (jugando con diferentes valor de min_samples), y a prtir de min_samples obtenemos una estimación del epsilon a utilizar.

# Epsilón adecuado
elbow_dbscan(X, 5)

# Cálculo
modelo_dbscan_moons = cluster.DBSCAN(eps=0.075, min_samples=5)
modelo_dbscan_moons.fit(X)

# Visualización
clusters_view(X, modelo_dbscan_moons.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

A diferencia de kmeans, DBSCAN al basarse en la densidad y no en la distancia en sí, conseguimos que en este caso, es decir, cuando los clústers no son convexos, identifique bien los grupos. Aún así vemos que se sigue habiendo outliers, hay pequeños puntos que no cumplen con el número de vecinos mínimos dentro del radio establecido.

En resumen, al ir calculado para cada punto si dentro del radio epsilon hay un número de vecinos mínimos, conseguimos que se identifiquen bien los clústers, ya que no se tiene en cuenta la distancia en sí de cada punto al centroide como sucede con kmeans, ya que en ese caso el centroide es difícil de establecer.

X, y = X_circles, y_circles

Implementación: prueba la implementación de DBScan jugando con los parámetros eps y min_samples para encontrar los grupos (y outliers) del dataset Circles.

Se ha utilizado un min_samples de 10 porque tras probar diferentes valores era lo que mejor funcionaba, y partir de ese valor se estima el valor de epsilon.

# Epsilón adecuado
elbow_dbscan(X, 10)

# Cálculo
modelo_dbscan_circles = cluster.DBSCAN(eps=0.1, min_samples=10)
modelo_dbscan_circles.fit(X)

# Visualización
clusters_view(X, modelo_dbscan_circles.labels_)

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

En este caso sucede lo mismo que en el ejercicio anterior, DBSCAN sí que es capaz de identificar bien los clústers ya que se basa en la densidad, es decir, no calcula la distancia entre un punto al centroide del clúster, sino que usa el radio (epsilon) y el número de vecinos mínimo (min_samples) para añadir un punto a un clúster.

Cabe destacar que en este caso también tenemos algunos outliers, pero son muy pocos. En resumen, DBSCAN consigue identificar de forma correcta los clústers.

1 c. Algoritmos jerárquicos

En este apartado se pide visualizar mediante un dendrograma la construcción progresiva de los grupos mediante un algoritmo jerárquico aglomerativo (estrategia bottom-up). Con ello se pretende encontrar un método gráfico para entender el comportamiento del algoritmo y encontrar los clusters deseados en cada dataset.

X, y = X_blobs, y_blobs

Implementación: prueba la implementación de clustering jerárquico de scipy probando distintos criterios de enlace o linkage permitiendo identificar los clusters subyacentes (mostrando su resultado) y su dendrograma para el dataset Blobs.
Puedes importar las librerías necesarias para ello.

import scipy.cluster.hierarchy as sch
from sklearn.cluster import AgglomerativeClustering

def dendrograma(X, method):
    plt.figure(figsize=(13, 8))
    plt.title("Dendrograma según " + method)
    sch.dendrogram(sch.linkage(X, method = method))
    plt.show() 

# Dendrograma
dendrograma(X, "single")

# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=2,
                                                    affinity='euclidean',
                                                    linkage='single')
modelo_aglomerativo_blobs.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)

# Dendrograma
dendrograma(X, "complete")

# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=4,
                                                    affinity='euclidean',
                                                    linkage='complete')
modelo_aglomerativo_blobs.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)

# Dendrograma
dendrograma(X, "average")

# Modelo
modelo_aglomerativo_blobs = AgglomerativeClustering(n_clusters=4,
                                                    affinity='euclidean',
                                                    linkage='average')
modelo_aglomerativo_blobs.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_blobs.labels_)

Análisis: Interpreta el dendrograma y comenta qué criterio de enlace se ha comportado mejor. ¿Por qué?

Para la realización de este ejercicio se ha calculado con tres tipos de enlace:

• Enlace simple: usa la distancia mínima entre elementos. Como se puede apreciar tanto en el dendrograma y los clústers que nos proporiona por salida no es el enlace adecuado, según el dendrograma hay 2 clústers y la ejecución del algoritmo vemos que hay dos clústers también, es decir, se produce el problema de cadena, el cual consiste en la unión de grupos solamente por tener elementos muy próximos entre los mismos.

• Enlace completo: usa la distancia máxima entre elementos. En este caso tanto en el dendrograma como en al ejecución del algoritmo nos proporciona los cuatros clústers que estamos buscando. Este tipo de enlace consigue eliminar el efecto cadena pero es sensible a los outliers. Aún así este enlace puede ser bueno para este problema.

• Enlace medio: es una mezcla entre el enlace simple y completo trantando de mitigar los problemas del efecto cadena y outliers, tanto en el dendrograma como en la ejecución del algiritmo nos indica de forma clara que lo mejor son cuatro clústers, por lo que este tipo de enlace también pueder ser útil para resolver el problema.

Como nota aclaratoria, el número de clústers ideales según el dendrograma viene determinado por la línea vertical más alta sin ser cortada horizontalmente, de tal forma que en ese línea hacemos un corte horizontal a todo el árbol y nos indica el número de clústers ideales para aplicar el algoritmo aglomerativo.

En resumen, tanto el enlace completo como el enlace medio pueden ser usados de forma idónea para resolver este problema.

X, y = X_moons, y_moons

Implementación: prueba la implementación de clustering jerárquico de scipy probando distintos criterios de enlace o linkage permitiendo identificar los clusters subyacentes (mostrando su resultado) y su dendrograma para el dataset Moons.
Puedes importar las librerías necesarias para ello.

# Dendrograma
dendrograma(X, "single")

# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
                                                    affinity='euclidean',
                                                    linkage='single')
modelo_aglomerativo_moons.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)

# Dendrograma
dendrograma(X, "complete")

# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
                                                    affinity='euclidean',
                                                    linkage='complete')
modelo_aglomerativo_moons.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)

# Dendrograma
dendrograma(X, "average")

# Modelo
modelo_aglomerativo_moons = AgglomerativeClustering(n_clusters=2,
                                                    affinity='euclidean',
                                                    linkage='average')
modelo_aglomerativo_moons.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_moons.labels_)

Análisis: Interpreta el dendrograma y comenta qué criterio de enlace se ha comportado mejor. ¿Por qué?

Para la realización de este ejercicio se han seguido usando los enlaces simple, completo y medio.

Respecto al enlace simple, vemos que en este caso el dendrograma nos indica de forma clara que dos clústers son lo ideal (porque la linea azul primera es la más grande y si partimos el diagrama horizontalmente por ahí nos quedan dos grupos). Al ejecutar el algoritmo con dos clústers vemos que realiza una perfecta clasificación, es decir, hay dos claros grupos diferenciados y al hacer uso de la distancia mínima entre elementos del grupo se consigue los clústers que estamos buscando, sin tener el efecto de cadena.

En cuanto al enlace completo, el dendrogrma nos indica que el número de clústers ideales pueden ser dos o cuatro, teniendo en cuenta la línea azul más a la izquierda tenemos 2 clústers, pero si nos fijamos en la línea verde más a la izquierda tendríamos 4 clústers. Como sabemos que lo ideal son dos, al ejecutar el algoritmo vemos que no consigue identificar bien los clústers, debido a que calcula la distancia máxima entre elementos.

Respecto al enlace medio, el dendrograma indica que el número ideal es 2 clústers y al ejecutar el algoritmo vemos que tampoco consigue obtener los clústers que buscamos.

En resumen, para este problema el enlace que mejor funciona es el enlace simple, ya que tiene en cuenta la distancia mínima entre elementos.

X, y = X_circles, y_circles

Implementación: prueba la implementación de clustering jerárquico de scipy probando distintos criterios de enlace o linkage permitiendo identificar los clusters subyacentes (mostrando su resultado) y su dendrograma para el dataset Circles.
Puedes importar las librerías necesarias para ello.

# Dendrograma
dendrograma(X, "single")

# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=2,
                                                    affinity='euclidean',
                                                    linkage='single')
modelo_aglomerativo_circles.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)

# Dendrograma
dendrograma(X, "complete")

# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=6,
                                                    affinity='euclidean',
                                                    linkage='complete')
modelo_aglomerativo_circles.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)

# Dendrograma
dendrograma(X, "average")

# Modelo
modelo_aglomerativo_circles = AgglomerativeClustering(n_clusters=5,
                                                    affinity='euclidean',
                                                    linkage='average')
modelo_aglomerativo_circles.fit_predict(X)

# Visualización
clusters_view(X, modelo_aglomerativo_circles.labels_)

Análisis: Interpreta el dendrograma y comenta qué criterio de enlace se ha comportado mejor. ¿Por qué?

Como era de esperar, este problema es similar al anterior por lo que el enlace que mejor funciona es el simple, básicamente porque tiene en cuenta la distancia mínima entre los elementos del grupo y como los clústers están bien diferenciados no se produce el efecto cadena. El dendograma de este enlace nos indica que el mejor número de clústers es 2, es decir, los clústers que estamos buscando.

En cuanto al enlace completo y enlace medio no son idóneos para este tipo de problemas, el primer enlace nos indica que el número de clústers tiene que ser 6 (hay que cortar por la penúltima línea vertical azul), y el segundo enlace nos indica que tiene que ser 5 grupos (hay que cortar por la tercera línea azul). Por lo tanto, los dendrogramas de estos enlaces nos están indicando que no son buenos enlaces para nuestro problema, ya que nos indican que lo ideal son 5 y 6 grupos cuando realmente son 2.

2. Aplicación: generación de imágenes con reducción de dimensionalidad (3 puntos)

Es posible aplicar una amplia variedad de algoritmos para la reducción de dimensionalidad. Para ello se empleará el dataset MNIST compuesto de miles de dígitos manuscritos del 0 al 9. Donde cada imagen se compone de 784 píxeles (imágenes de 28 x 28), por lo que se parte de un número alto de dimensiones.

X, y = sklearn.datasets.fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False)

Por lo que cada muestra (las 70k filas del dataset) se componen de 784 dimensiones:

X.shape

(70000, 784)

fig, axis = plt.subplots(1, 10, figsize=(12, 6))
for i, ax in enumerate(axis):
    ax.imshow(X[y == str(i)][0].reshape(28, 28), cmap='gray_r')
    ax.set_title(str(i))
    ax.axis('off')
plt.tight_layout()

Si cada algoritmo obtiene resultados distintos a la hora de reducir la dimensionalidad, ¿qué representación es más fiel a la distribución original?

Antes de reducir las 784 dimensiones originales de cada muestra a 2 para poder visualizarlas en 2 dimensiones, es muy útil conocer, o al menos intuir, la estructura en alta dimensionalidad de los datos.

Para ello se puede hacer uso del dendrograma como heurística para conocer la disposición original de los datos y comprobar si la proyección es similar a lo mostrado por el dendrograma.

Implementación: aprender una proyección a 2 dimensiones de las muestras de X con PCA y proyectar el conjunto X a dos dimensiones. Después visualizarlo en un scatter plot. Utiliza las etiquetas de y (el número manuscrito al que se corresponde cada muestra), en el parámetro label (en la llamada a scatter) y la función legend en la visualización para saber la clase correspondiente a cada punto e interpretar el resultado de la reducción de dimensionalidad y poder interpretar el resultado de la proyección.

# Creación del PCA
modelo_pca = decomposition.PCA(n_components=2)
X_pca = modelo_pca.fit_transform(X)

# Visualización
# https://www.kaggle.com/parulpandey/part1-visualizing-kannada-mnist-with-pca
y_numerico = [int(i) for i in y] # Convertimos a número las etiquetas

plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(X_pca[:, 0], X_pca[:, 1], s=5, c=y_numerico, cmap='Spectral')
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('Análisis de componentes principales', fontsize=14)
plt.xlabel('Primer componente principal', fontsize=14)
plt.ylabel('Segundo componente principal', fontsize=14)

Text(0, 0.5, 'Segundo componente principal')

Análisis: ¿Qué puedes interpretar de la proyección? ¿Las clases han quedado visiblemente separadas? ¿Por qué?

Lo primero de todo, las clases no han quedado visiblemente separadas del todo, es decir, sí que se pueden apreciar patrones gracias a los colores, pero están todos los puntos juntos y es difícil de diferenciar en la parte central de la representación.

Uno de los motivos por los que no queda del todo bien diferenciadas las clases puede ser porque estamos reduciendo las dimensiones originales (784) a dos dimensiones, es decir, el modelo no es capaz de representar toda la variabilidad del mismo ya que hay mucha diferencia entre las dimensiones originales a dos dimensiones (se pierde información). Ésto se podría solventar haciendo un estudio de cuántas dimensiones son ideales para este problema, y a partir de ahí generar el modelo con dichas dimensiones. Otro punto pueder por que haya demasiada información en el dataset, es decir, demasiados registros.

Aunque se ha comentado que es difícil de interpretar en profundiad la representación, sí que podemos ver patrones en los datos. Por ejemplo, vemos que el número 0 (rojo oscuro) y el número 7 (verde azulado) están opuestos el uno del otro, esto tiene toda la lógica del mundo y es que los píxeles usados para representar ambos son muy diferentes. Otro caso sería el del 1 (rojo) junto con el 0 y el 7, es decir, el uno usa píxeles muy parecidos al número 7 por lo que ambos están más juntos en la representación que respecto al número 0. También vemos que en el medio de la representación hay muchos números que usan píxeles parecidos, como el 4, el 5, el 6.

En resumen, aunque las clases no quedan del todo definidas sí que lo hacen de forma suficiente como para poder detectar patrones en los datos.

En la gráfica anterior cada punto representa una muestra en 2 dimensiones. Con PCA es posible invertir la transformación para que, a partir de cada punto 2d, se obtenga de nuevo (aproximadamente) la imagen original (784 dimensiones).

Por lo que es posible "generar" nuevas imágenes eligiendo puntos aleatoriamente del plano 2d, y pedirle al modelo PCA aprendido que invierta la transformación para obtener las "teóricas" imágenes que habrían sido proyectadas a esos puntos del espacio proyectado.

Implementación: calcular el máximo y mínimo de cada una de las dos dimensiones y, para cada una de ellas, generar una secuencia de 10 valores con igual separación.

# Mínimos y máximos
c1_min = min(X_pca[:,0])
c1_max = max(X_pca[:,0])
c2_min = min(X_pca[:,1])
c2_max = max(X_pca[:,1])

# Secuencias
c1_sec = np.linspace(start=c1_min, stop=c1_max, num=10)
c2_sec = np.linspace(start=c2_min, stop=c2_max, num=10)

Con las dos secuencias de 10 (una por cada dimensión del plano de proyección) valores es posible combinar los puntos de ambas secuencias para generar 100 combinaciones (puntos 2d) que teselan el plano sobre el que PCA ha proyectado las muestras.

Implementación: invertir la transformación para cada uno de los 100 puntos y visualizar su imagen asociada en una matriz de 10 x 10 imágenes (tratando de preservar su posición en el espacio proyectado).

# Método que devuelve las secuencias ordendas en el espacio proyectado
def sorted_sec(c1_sec, c2_sec):
    secuencia = []
    c2_sec_aux = c2_sec
    
    while len(c2_sec_aux) > 0:
        maximo = max(c2_sec_aux)
        pos_max = np.argmax(c2_sec_aux)
        
        for c1 in c1_sec:
            secuencia.append([c1,maximo])
            
        c2_sec_aux = np.delete(c2_sec_aux, pos_max)
        
    return np.array(secuencia)

secuencia = sorted_sec(c1_sec, c2_sec)

fig, axis = plt.subplots(10, 10, figsize=(10, 10))
row = 0
col = 0
for i in range(len(secuencia)):
    if i % len(c1_sec) == 0 and i > len(c1_sec) - 1: # Si empieza una nueva fila y no es la primera
        col = 0
        row += 1
    number = modelo_pca.inverse_transform(secuencia[i])
    axis[row,col].imshow(number.reshape(28, 28), cmap='gray_r')
    axis[row,col].axis('off')
    col += 1
plt.tight_layout()

Análisis: ¿Qué puedes interpretar de las imágenes reconstruidas / interpoladas? ¿Genera números o transiciones entre números visualmente creíbles? ¿Por qué?

Sí que se puede interpretar las imágenes reconstruidas, ya que éstas se corresponden con el espacio proyectado anteriormente, pero no son del todo claras dichas representaciones.

Por ejemplo, en la parte superior vemos que el número que más se intuye es el 9 y a medida que nos desplazamos a hacia la derecha ese nueve se va transformado hasta llegar a una especie de 0 sin ser el 0. Si partimos de la zona media de la representación vemos que el número que más se intuye es el 1 y éste se va transformando por 3/5/8 hasta llegar al 0. Y si nos fijamos en la zona baja de la representación vemos como una especie de 1 sin llegar a ser uno, se va transformando por 2 y 3 hasta llegar a una especie de 2 y 0.

Es decir, vemos que sí se generan números y transiciones creíbles según el espacio proyectado que se ha graficado al principio del ejercicio, pero una cosa es creíble y otra entendible, ya que de lo segundo se requiere de una gran abstracción. Esto sucede porque hemos generado diferentes puntos dentro del espacio proyectado, para ser más exactos 100, los cuales de forma equidistante nos permiten a partir de PCA obtener su representación, lógicamente estos puntos no son del dataset original, por lo que pueden quedar bien representados o ser una mera transición entre números.

En resumen, vemos que sí que se consigue obtener números o transiciones no muy bien definidas entre números a partir de puntos equidistantes en el plano gracias a poder calcular la transformación inversa.

Análisis: ¿Podría haberse conseguido lo mismo con t-SNE? ¿Por qué?

No, no se podría haber conseguido lo mismo con t-SNE, esto se debe a que nos permitiría realizar la primera parte del ejercicio pero no la segunda, es decir, t-SNE nos permitiría visualizar las diferentes clases en el espacio proyectado, pero no nos permitiría calcular la transformación inversa, por lo tanto no conseguiríamos que a partir de unos puntos equidistantes del plano se pudieran obtener los números o las transiciones entre los mismos.

Implementación: aprender una proyección a 2 dimensiones de las muestras de X con UMAP (con los parámetros por defecto) usando la librería umap-learn y proyectar el conjunto X a dos dimensiones. Después visualizarlo en un scatter plot. Utiliza las etiquetas de y (el número manuscrito al que se corresponde cada muestra), en el parámetro label (en la llamada a scatter) y la función legend en la visualización para saber la clase correspondiente a cada punto e interpretar el resultado de la reducción de dimensionalidad y poder interpretar el resultado de la proyección.

# Creación del UMAP
modelo_umap = umap.UMAP(n_components=2, random_state=0)
X_umap = modelo_umap.fit_transform(X)

# Visualización
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(X_umap[:, 0], X_umap[:, 1], s=5, c=y_numerico, cmap='Spectral')
plt.gca().set_aspect('equal', 'datalim')
plt.colorbar(boundaries=np.arange(11)-0.5).set_ticks(np.arange(10))
plt.title('Análisis de UMAP', fontsize=14)
plt.xlabel('Primer componente principal', fontsize=14)
plt.ylabel('Segundo componente principal', fontsize=14)

Text(0, 0.5, 'Segundo componente principal')

Análisis: ¿Qué puedes interpretar de la proyección? ¿Las clases han quedado visiblemente separadas? ¿Por qué?

En este caso UMAP a nivel visual se comporta mejor que PCA, ya que sí que quedan visiblemente separadas cada una de las clases.

Vemos que podemos seguir identificando patrones, ya que hay números que utilizan unos píxeles similares como el grupo formado por los números 7, 9 y 4. Otro grupo formado por los números 8, 5 y 3, y el resto de números que son grupos independientes, aunque por ejemplo el número 6 está más cerca del 0 que del 1.

La principal diferencia entre PCA y UMAP es qu el primero es un algoritmo lineal y el segundo no, es decir, PCA está basado en la factorización de matrices mientras que UMAP está basado en gráfos, es por ello que UMAP consigue identificar mejor a nivel visual las clases.

Implementación: al igual que anteriormente con PCA, calcula el máximo y mínimo para cada una de las dos dimensiones e invierte la transformación con el modelo aprendido por UMAP para cada uno de los 100 puntos y visualizar su imagen asociada en una matriz de 10 x 10 imágenes (tratando de preservar su posición en el espacio proyectado). Consejo: la inversión de la transformación en UMAP es más costosa computacionalmente que con PCA, por lo que recomiendo que sólo la invoques una vez con las 100 muestras en lugar de hacer 100 llamadas (una por cada muestra). Esto reduce drácticamente el tiempo de ejecución.

# Mínimos y máximos
c1_min = min(X_umap[:,0])
c1_max = max(X_umap[:,0])
c2_min = min(X_umap[:,1])
c2_max = max(X_umap[:,1])

# Secuencias
c1_sec = np.linspace(start=c1_min, stop=c1_max, num=10)
c2_sec = np.linspace(start=c2_min, stop=c2_max, num=10)

# Secuencia ordenada para respetar el espacio proyectado
secuencia = sorted_sec(c1_sec, c2_sec)

# Números de la representación
representaciones = modelo_umap.inverse_transform(secuencia)

# Representación
fig, axis = plt.subplots(10, 10, figsize=(10, 10))
row = 0
col = 0
for i in range(len(representaciones)):
    if i % len(c1_sec) == 0 and i > len(c1_sec) - 1: # Si empieza una nueva fila y no es la primera
        col = 0
        row += 1
    number = representaciones[i]
    axis[row,col].imshow(number.reshape(28, 28), cmap='gray_r')
    axis[row,col].axis('off')
    col += 1
plt.tight_layout()

Análisis: ¿Qué puedes interpretar de las imágenes reconstruidas / interpoladas? ¿Genera números o transiciones entre números visualmente creíbles? ¿Por qué?

En este caso vemos que las imágenes reconstruidas son perfectamente entendibles, más que lo que sucedía con PCA, es decir, nos genera principalmente números aunque también hay alguna que otra transición.

EL porqué es sencillo, UMAP trata de "clasificar" los datos además de reducir la dimensionalidad, con ello me refiero a que busca que los grupos de datos estén lo más diferenciados posible, mientras que PCA busca proyectar cada punto en el hiperplano buscando la máxima variabilidad independientemente de si los grupos están más cerca o no. Esto hace que cuando hagamos la transoformación inversa UMAP funcione mejor, ya que hay una mayor separación/diferencia entre los grupos, mientras que en PCA creaba un grupo muy grande en el que estaban todos los datos.

3. Aplicación: identificación de puntos de interés turísticos (3 puntos)

En este ejercicio se busca automatizar la localización de lugares turísticos a través de los metadatos de las fotografías de flickr.

Para ello se provee junto a la PEC el dataset: barcelona.csv. Ya que se pide encontrar los puntos de mayor interés turístico de esta ciudad.

Opcional: si quieres hacerlo para otra región

Pero si quieres hacerlo para otra parte del mundo, puedes descargarte el dataset completo aquí y descomprime para extraer el CSV.

Para seleccionar las coordenadas de la zona de interés puedes usar la opción Export manual de OpenStreetMaps.

Por último, para filtrar los datos que se corresponden a la zona deseada puedes usar el programa AWK mediante la siguiente línea:

awk -F"," 'NR == 1 {print $5","$6} (NR > 1 && $5 > 41.3560 && $5 < 41.4267 && $6 > 2.1300 && $6 < 2.2319) {print $5","$6}' photo_metadata.csv

$5 hace referencia a la latitud, y $6 a la longitud. Sustituye los valores mínimo y máximo para obtener los datos de localización referentes a tu área de interés.

geo_df = pd.read_csv('./data/barcelona.csv', header=0)
geo_df.sample(5)

	latitude	longitude
12408	41.404778	2.174541
2740	41.391833	2.165167
1701	41.368942	2.147912
11082	41.403145	2.212886
4379	41.386500	2.155999

Implementación: siempre que tratamos un problema real, es necesario entender los datos a tratar. Visualiza las localizaciones de las fotografías mediante un scatter plot. Prueba distintos parámetros de tamaño (size) s, y opacidad alpha hasta conseguir un resultado fácil de interpretar.

def show_barcelona(df):
    plt.figure(figsize=(13,8), facecolor="white")
    plt.scatter(df.latitude, df.longitude, s=10, alpha=.25)
    plt.title("Lugares turísticos de Barcelona")
    plt.xlabel("Latitud")
    plt.ylabel("Longitud")

show_barcelona(geo_df)

Análisis: tras haber probado los algoritmos de agrupamiento en el ejercicio 1. ¿Qué algoritmo crees que sería más adecuado tras visualizar los datos? ¿Por qué?

El algoritmo más adecuado para aplicar clustering sobre el problema es DBSCAN, es decir, un algoritmo de densidad.

Es verdad que se podría hacer uso de otros algoritmos como por ejemplo KMEANS, pero éste se basa en el cálculo de centroides y asigna cada punto al centroide más cercano, esto puede ser un inconveniente porque dado este problema podría considerar varios clústers como uno mismo. Además, tanto KMEANS como los algortimos jerárquicos necesitan pasarle por parámetro al algoritmo el número de clústers, cosa que en este ejercicio se desconoce, y lo que se busca es conocer eso número de clústers, es decir, cuántos puntos turísticos hay en Barcelona. También ha que destacar que DBSCAN es el único en detectar outliers, siendo éstos los puntos menos turísticos.

En resumen, el algoritmo DBSCAN es el que mejor funciona según nuestro problema ya que se basa en la estrucutra/densidad de los datos, de esta forma podemos conseguir qué clústers hay en nuestro problema y cómo son.

Implementación: para prototipar el modelado primero se recomienda elegir un subconjunto de los datos que sea representativo. Selecciona una muestra del DataFrame original y visualiza como en el punto anterior para comprobar su similitud.

geo_df_muestra = geo_df.sample(frac=.3, random_state=0)
show_barcelona(geo_df_muestra)

Implementación: ajusta el algoritmo de clustering elegido para encontrar los distintos grupos sobre el conjunto reducido, y visualiza el resultado coloreando cada punto en base al grupo al que pertenece. Como pista, alrededor de 20 clusters es un número razonable, y es posible darles un color distinto a cada uno con el colormap: tab20.

# Estmiación de epsilon
elbow_dbscan(geo_df_muestra, 20)

# Modelo
modelo_dbscan_barcelona = cluster.DBSCAN(eps=0.0025, min_samples=20)
modelo_dbscan_barcelona.fit(geo_df_muestra)

# Visualizamos
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(geo_df_muestra.latitude, geo_df_muestra.longitude, s=10, alpha=.25, c=modelo_dbscan_barcelona.labels_, cmap="tab20")
plt.title("Lugares turísticos de Barcelona")
plt.xlabel("Latitud")
plt.ylabel("Longitud")

Text(0, 0.5, 'Longitud')

Implementación: si has usado un método de clustering que permite la detección de outliers. Representa sólo los puntos que no ha considerado outliers, es decir, los que pertenecen a algún cluster.

# Limpieza, si el clúster es -1 entonces es un outiler
geo_df_muestra["Target"] = modelo_dbscan_barcelona.labels_
geo_df_clean = geo_df_muestra[geo_df_muestra["Target"] != -1]

# Visualizamos
plt.figure(figsize=(13,8), facecolor="white")
plt.scatter(geo_df_clean.latitude, 
            geo_df_clean.longitude, 
            s=10, 
            alpha=.25, 
            c=geo_df_clean.Target, 
            cmap="tab20")
plt.title("Lugares turísticos de Barcelona")
plt.xlabel("Latitud")
plt.ylabel("Longitud")

Text(0, 0.5, 'Longitud')

Análisis: interpreta cual es el lugar que representa cada cluster (si encuentras una asociación lógica).

Antes de hacer el análisis vamos a calcular el punto medio de cada clúster, el cual nos va a servir tanto para este ejercicio como para el siguiente:

# Método que devuelve la media de la longitud y latitud para cada clúster
def get_meanCluster(df):
    dicMeanCluster = {}
    
    for i in np.unique(df.Target):
        latitud = np.mean(geo_df_clean[geo_df_clean.Target == i].latitude)
        longitud = np.mean(geo_df_clean[geo_df_clean.Target == i].longitude)
        dicMeanCluster[i] = [latitud, longitud]
        
    return dicMeanCluster

# Método que imprime la media de la longitud y latitud de cada clúster
def print_MeanCluster(dic):
    print("\t\t\tLatitud\t\tLongitud")
    for k,v in dic.items():
        print("Clúster " +
              str(k) +
              "\t->\t" +
              str(round(v[0], 7)) +
              "\t" +
              str(round(v[1], 7)))      

mean_clusters = get_meanCluster(geo_df_clean)
print_MeanCluster(mean_clusters)

			Latitud		Longitud
Clúster 0	->	41.3851449	2.1749851
Clúster 1	->	41.3713111	2.150776
Clúster 2	->	41.3970314	2.1904491
Clúster 3	->	41.4132301	2.1991785
Clúster 4	->	41.4015635	2.1682176
Clúster 5	->	41.4062836	2.2015134
Clúster 6	->	41.3714757	2.1402785
Clúster 7	->	41.4036251	2.1900497
Clúster 8	->	41.4107145	2.2215941
Clúster 9	->	41.4038759	2.1743124
Clúster 10	->	41.4144028	2.1526811
Clúster 11	->	41.3645843	2.1642466
Clúster 12	->	41.3753117	2.1570397
Clúster 13	->	41.3805483	2.1546133
Clúster 14	->	41.3970029	2.1946244
Clúster 15	->	41.3984371	2.2124015
Clúster 16	->	41.3983408	2.1824438
Clúster 17	->	41.3783561	2.1613572
Clúster 18	->	41.3688051	2.1608188
Clúster 19	->	41.4116274	2.1743459

Para interpretar cada clúster se ha calculado el punto medio de cada uno y se ha buscado dicha latitud y longitud en Google (por desconocimiento de la ciudad de Barcelona). Es verdad que para muchos de los puntos calculados sí que aportan sitios significativos, pero para otros no tanto, seguramente haya más que puntos zonas/barrios de interés.

A partir de la anterior visualización y la búsqueda en Google hemos obtenido esta información:

• Clúster 0: Azul oscuro es el centro histórico de Barcelona (Catedral de Barcelona, barrio Gótico).
• Clúster 1: Azul violeta es la zona de Montjuic.
• Clúster 2: Naranja es el el barrio Poblenou, la zona del teatro nacional.
• Clúster 3: El punto naranja solitario es díficil de interpretar, porque no hay ningún museo, catedral, nada, y hace referencia la Gran vía de las Cortes Catalanas en el barrio de San Martín.
• Clúster 4: El verde oscuro son los alrededores de la sagrada familia.
• Clúster 5: El verde claro es el parque de Poblenou.
• Clúster 6: El rojo oscuro son los alrededores de Montjuic.
• Clúster 7: El rojo claro es la torre Agbar.
• Clúster 8: El morado oscuro la zona del Auditori Forum, museo de ciencias.
• Clúster 9: El morado clarito es la sagrada familia.
• Clúster 10: El marrón oscuro es el Park Güel.
• Clúster 11: El marrón claro el castillo de Montjuic.
• Clúster 12: Azul violeta es la plaza España.
• Clúster 13: El rosa los jardínes de Emma.
• Clúster 14: Gris oscuro representa la zona del parque Mercé Plantada.
• Clúster 15: El gris claro representa la playa de la Nova Mar Bella.
• Clúster 16: El verde pistacho oscuro hace referencia a la zona de la Monumental.
• Clúster 17: El verde pistacho claro es la zona del mercado de Sant Antoni.
• Clúster 18: El azul claro representa la zona de la fundación Joan Miró.
• Clúster 19: El azul claro al lado del morado claor hace referencia al Recinte Modernista de Sant Pau.

Cabe destacar que la visualización es un poco costosa de entender, esto se debe a que la parte más alta hace referencia al mar, esta parte se puede intuir así, pero si al "norte" está el mar mirando Google Montjuic tendría que estar a la derecha, y vemos que esto no es así, es decir, Montjuic está a la izquierda. Esto nos da a entender que la visualización no es la realidad sino que es un modo espejo de la realidad.

OPCIONAL Implementación: representa los puntos sin ruido sobre un mapa utilizando la librería Smopy. Para facilitar la interpretación, puedes representar cada cluster como el punto medio de todos los puntos que lo conforman.

import smopy

mapa = smopy.Map((41.36, 2.13, 41.42, 2.23), z=12)

lista_x = []
lista_y = []
for k, v in mean_clusters.items():
    x, y = mapa.to_pixels(v[0], v[1])
    lista_x.append(x)
    lista_y.append(y)
ax = mapa.show_mpl(figsize=(15,15))
ax.scatter(lista_x, lista_y, s=350, c="purple")

<matplotlib.collections.PathCollection at 0x25e21dc97c0>